Coarsening of two dimensional XY model with Hamiltonian dynamics: 
Logarithmically divergent vortex mobility 
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We investigate the coarsening kinetics of an XY model defined on a square lattice when the 
underlying dynamics is governed by energy-conserving Hamiltonian equation of motion. We find 
that the apparent super-diffusive growth of the length scale can be interpreted as the vortex mobility 
diverging logarithmically in the size of the vortex-antivortex pair, where the time dependence of the 
characteristic length scale can be fitted as L(t) ~ ((t + to) ln(t + to)) 1 / 2 with a finite offset time to- 
This interpretation is based on a simple phenomenological model of vortex-antivortex annihilation 
to explain the growth of the coarsening length scale L(t). The nonequilibrium spin autocorrelation 
function A(t) and the growing length scale L(t) are related by A(t) ~ L~ x (t) with a distinctive 
exponent of A ~ 2.21 (for E = 0.4) possibly reflecting the strong effect of propagating spin wave 
modes. We also investigate the nonequilibrium relaxation (NER) of the system under sudden heating 
of the system from a perfectly ordered state to the regime of quasi-long-range order, which provides a 
very accurate estimation of the equilibrium correlation exponent rj(E) for a given energy E. We find 
that both the equal-time spatial correlation C„ r (r,t) and the NER autocorrelation A nr (t) exhibit 
scaling features consistent with the dynamic exponent of z nr = 1. 

PACS numbers: 64.60. Ht,64.60.Cn,75. 10. Hk,75.40.Gb 



I. INTRODUCTION 

Phase ordering kinetics refers to a study of nonequi- 
librium dynamic processes of statistical mechanical sys- 
tems approaching toward the thermal equilibrium of or- 
dered states with broken symmetry. The system is usu- 
ally quenched from a random disordered state to a low- 
temperature ordered state below the transition tempera- 
ture 0, Q • The ever-nonequilibrium approach to ordered 
state, in the thermodynamic limit, is usually character- 
ized by the dynamic scaling accompanied by power-law 
growth of typical length scales of the system. These 
power law exponents depend on the spatial dimension 
and the symmetry of the ground states. More impor- 
tantly these power laws also depend on the nature of the 
dynamic protocols such as whether the order parameter 
is conserved or not during the phase ordering processes. 
Usually, characteristic topological defects, such as point 
vortices or domain walls, are generated in the initial dis- 
ordered state, depending on the dimension of the order 
parameter, and the phase ordering processes can also be 
interpreted as the process of annihilation and decay of 
these topological defects: the motion of these defects 
(and their annihilation) determines the characteristics of 
the coarsening processes. 

Most of the computational works on phase ordering 
kinetics have been based on a purely dissipative dynam- 
ics. In this case, we expect that the motion of the vor- 
tices would be also purely dissipative. Energy scaling 
arguments in combination with conservation laws could 
determine the growth laws for most of the 0(n) spin sys- 
tems in d dimensions [2| (except for the special case of 
d = n = 2). 



In reality, however, there exist various systems exhibit- 
ing dynamic processes that cannot be described solely by 
dissipative dynamics. For example, inertial effects would 
be present in the case of magnetic spin systems, where 
the spins are influenced by neighboring spins via preces- 
sion interaction terms that are energy conserving 
that are often called reversible mode coupling. Actually, 
among the eight model systems classified by Hohenberg 
and Halperin 0, @] for describing the dynamic critical 
phenomena, five of them (models E, F, G, H, and J) in- 
clude reversible mode couplings. 

We may also take a microcanonical approach to the 
nonequilirium dynamics by modeling the dynamics of the 
system based on a purely Hamiltonian dynamics with en- 
ergy conservation. In relation to this, there have been 
some research efforts in recent years toward investigating 
the statistical mechanical behavior of lattice spin systems 
by solving directly the Hamiltonian equation of motion 
associated with the spin system [8j]. For example, equi- 
librium phase transitions have been found to be related 
to a topological change in the dynamics ■ In terms 

of nonequilibrium phase ordering dynamics, Zheng inves- 
tigated via Hamiltonian dynamics the coarsening of the 
Ising-type <fr model on a two-dimensional square lattice 
with nonconserved order parameter [TH [l6| . 

In this work, we study numerically the coarsening pro- 
cess of an XY model on a square lattice using Hamil- 
tonian dynamics. This work is built upon a previous 
preliminary work [l7| on the same model system, which 
we think is incomplete, especially in terms of the analysis 
on the growth law. 

In equilibrium, the XY model exhibits a Berezinskii- 
Kosterlitz-Thouless (BKT) transition at Tbkt due to the 
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unbinding of vortex-antivortex pairs 18]. Below Tbkt, 
the system has a quasi-ordered phase that is character- 
ized by a power-law decay of the order parameter cor- 
relation function for long distances. The critical expo- 
nent governing the power-law decay decreases continu- 
ously down to zero temperature: the system is critical at 
equilibrium for all non-zero temperatures below Tbkt- 

The phase ordering dynamics of the XY model has 
been studied for quite a long time by many groups [l9l - 
l29l ]. It is now agreed that, in the phase ordering dy- 
namics of the XY model with non-conserving order pa- 
rameter (via Monte Carlo simulations or Langevin dy- 
namics methods), the growing l eng th scale £mc(^) ex ~ 
hibits a logarithmic correction 30] to diffusive growth 

1 /2 

as Luc(t) ~ (t/logf) . Since the coarsening in 
XY model is dominated by the annihilation of vortex- 
antivortex pairs generated in the random initial states, 
we may interpret the growth law in terms of the mo- 
bility (or friction) characteristics of the vortices. Here, 
the logarithmic correction can be attributed to a log- 
arithmic divergence (in the system size) of the effective 
friction constant of a moving vortex in the dissipative dy- 
namics with a non-conserving order parameter. Another 
system that can be described by XY model at equilib- 
rium is an array of superconducting Josephson junctions. 
In the cases of superconducting-normal-superconducting 
(SNS) junction arrays, the dynamics is often described 
by resistively-shunted junction (RS J) model. This model 
is characterized by a Laplacian type of coupling between 
time derivatives of superconducting phases of neighbor- 
ing islands. Due to this particular type of dissipative 
coupling in Josephson junction arrays, there is no loga- 
rithmic correction in the growing length scale that cor- 
responds to a finite friction constant of a moving vortex 
31]. 

One further experimental system that can be described 
by XY model is a thin film of supcrfluid helium. The 
dynamics of this system is characterized by a reversible 
mode-coupling term. In recent simulations on the coars- 
ening dynamics of XY model with reversible mode- 
coupling (MCXY), it was shown that the growth law 
exhibits a positive logarithmic correction which is oppo- 
site to the coarsening of ordinary XY model with purely 
dissipative dynamics and non-conserving order parame- 
ter [5j. This system apparently exhibits a growth law 
with an exponent that is a little larger than the diffusive 
growth exponent of 1/2. But this was re-analyzed and 
fitted in terms of a (positive) logarithmic correction to a 
diffusive growth. The growth law could also be under- 
stood in terms of a simple vortex-antivortex annihilation 
model. These features reminded us of the growth law 
with apparent super-diffusive exponents in the coarsen- 
ing of XY model with Hamiltonian dynamics (l7j . 

Now, in this work, we perform more extensive simula- 
tions on the Hamiltonian coarsening of XY model and 
carry out an analysis of the growth law in more details 
in order to determine whether these similarities between 
MCXY and Hamiltonian XY models are indeed valid. 



Due to the Hamiltonian dynamics nature of our system, 
the total energy is conserved, and we specify the initial 
states by their total energy divided by the system size, 
i.e., the per-site energy E. These initial states with spec- 
ified energies are prepared with a special Monte Carlo 
algorithm [l7[ ■ The initial rotational velocities of the ro- 
tors are taken to be zeros, ensuring that the initial con- 
figuration is maximally disordered within the constraint 
of the fixed energy. 

Since we begin with zero kinetic energy (zero rotational 
velocities for all rotors), the Hamiltonian dynamics of 
the system will generate kinetic energy taken from the 
potential energy. As time proceeds, we can expect that 
(for a per-site energy that is low enough to correspond to 
a low-temperature quasi-ordered phase), the system will 
evolve toward some equilibrium stationary state which 
can be considered as corresponding to the ordinary ther- 
mal equilibrium state with critical quasi-long-range or- 
der. In other words, the kinetic energy that is generated 
during the course of coarsening process acts as a thermal 
bath for the system. We investigate the time-dependent 
spin configuration for the system in terms of spatial or- 
dering and relaxation of the vortex numbers, etc., in anal- 
ogy to the conventional dissipative coarsening systems. 

We find that the equal-time spatial correlation func- 
tions satisfy critical dynamic scaling 

C( r> t)=r-^g(r/L(t)), (1) 

with a spatial correlation exponent tj(E) (which is in- 
creasing in the per-site energy E) and a growing length 
scale L(t) that grows typically as L(t) ~ t 1 '", where z 
is the dynamic exponent. In the late time region, we 
find that the length scale L(t) grows with an exponent 
1/z that is apparently a little larger than the diffusive 
exponent 1/2. However, from more detailed analysis we 
could confirm that the growth of the length scale can 
be interpreted as the vortex mobility diverging logarith- 
mically in the size of the vortex-antivortex pair. The 
growth of characteristic length scale can be fitted as 
L(t) ~ ((t + t )\n(t + h)) 1/2 with a finite offset time 
to. This offset time to corresponds to a finite length scale 
Lo at t = due to the finite average spacing between 
vortices and antivortices at the initial time. 

We could also compare the nonequilibrium spin auto- 
correlation function A{t) of the order parameter in the 
coarsening process and the growing length scale L(t). 
Here, we found that A(t) and the L(t) are related by 
a scaling exponent A with A(t) ~ L~ x (t) where A ~ 2.21 
(at E = 0.4) with only very weak monotonic increase 
on the value of E. This is rather close to that in the 
case of MCXY or (soft-spin) model E @. While, on 
the other hand, this value of A exponent is quite differ- 
ent from that of usual dissipative coarsening dynamics 
of XY model (or 0(2) model) with nonconserved order 
parameter where the exponent of A ~ 1.17 was predicted 
theoretically and also observed in numerical simulations 
[2R I32l435j . This implies that the system loses its mem- 
ory much more quickly in the Hamiltonian dynamics or in 
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the case of MCXY than in the case of simple dissipative 
dynamics. 

Our simulation results indicate that the phase order- 
ing kinetics of Hamiltonian XY model belongs to the 
same universality class as that of the Langevin dynam- 
ics of XY model with a reversible mode-coupling and 
also that of (soft-spin) model E dynamics. This seems 
to be due to the existence of the same symmetry and 
conservation law in the two systems. That is, the total 
rotational momentum is conserved in the Hamiltonian 
dynamics which corresponds to the conservation of the 
third component (or the conjugate momentum) in the 
model E dynamics. Due to the conservation of the rota- 
tional momentum, propagating spin wave modes appear 
[36l [37| . This propagating spin wave modes might inter- 
act with the vortices and antivortices such that the mo- 
bility diverges logarithmically as the distance between 
vortices and antivortices increases. As in the previous 
work on the phase ordering dynamics of MCXY model, 
here also we attempt to fit the growth of the length scale 
with a simple phenomenological dynamic model of coars- 
ening based on vortex-antivortex annihilation. 

In order to analyze the critical dynamic scaling of the 
equal-time spatial correlation of the order parameter, we 
need an accurate estimate of the spatial correlation ex- 
ponent i](E). For accurate measurement of r)(E), we em- 
ploy an analog of the so-called nonequilibrium relaxation 
(NER) method. Here one starts with a perfectly ordered 
initial state but with random finite rotational velocities 
such that the total kinetic energy per site is equal to 
E. Beginning with this initial state, the system evolves 
toward a steady dynamic state corresponding to a criti- 
cal equilibrium state. By analyzing the approach to the 
steady state in terms of the equal-time spatial correlation 
function of the order parameter, we can estimate more 
accurately the equilibrium correlation exponent rj(E). In 
this case of NER relaxation, we can also obtain approx- 
imate analytic expression for the equal-time spatial cor- 
relation and autocorrelation function by using spin-wave 
approximation (see APPENDIX A), which agrees quali- 
tatively well with simulation results. 



II. THE HAMILTONIAN XY MODEL AND 
SIMULATION METHODS 



The dynamic equation of the hard-spin XY model on 
a square lattice (of linear size N) can be obtained from 
the following Hamiltonian given by 



J 



<ij> 



angular momentum associated with the angular variable 
9i of planar spin at site i. This model appears naturally 
in easy-plane ferromagnets or superfluid helium with the 
conserved variable m corresponding to the z-component 
of the spin or the density of the superfluid in the limit 
of negligible thermal noise (e.g., negligible phonon effect) 



Previous works [11|, UM have shown that a BKT tran- 
sition occurs at the value of the per-site energy around 
Ebkt — 1-0. In this work, we perform simulations 
of Newtonian dynamics based on the above Hamilto- 
nian with initial states of zero kinetic energy and a pre- 
specified potential energy E (but otherwise with random 
initial phases) . Specified potential energy values are cho- 
sen from the region that corresponds to below the BKT 
transition energy {Ebkt — 1-0). 

Unlike the case of dissipative systems with Langevin 
dynamics or Monte Carlo dynamics where dissipative re- 
laxation is incorporated in a natural way, here in our case 
of Hamiltonian dynamics, the total energy (sum of ki- 
netic and potential energy) is conserved. However, since 
we begin the Hamiltonian dynamics with zero kinetic en- 
ergy with a given specified maximal potential energy, we 
expect that the system will evolve in such a way as to de- 
crease (on the average) the potential energy and increase 
the kinetic energy. Eventually, in the long-time limit, the 
system will reach a certain kind of near-equilibrium dy- 
namic state where the average potential energy and the 
kinetic energy do not vary appreciably. Of course, rigor- 
ously speaking, in the thermodynamic limit, the system 
will never reach a complete equilibrium state in a finite 
time period. 

Therefore, we can consider the kinetic part of the sys- 
tem acting as a kind of thermal heat bath that absorbs 
the "dissipated heat" coming from the decreasing poten- 
tial energy. We can, thus, investigate the time-dependent 
spin configurations for the system in terms of spatial or- 
dering and relaxation of vortex numbers, etc. We gener- 
ate the initial states with specified potential energies in 
the following manner. To begin with, we take a random 
spin configuration. Then we employ a Monte Carlo an- 
nealing algorithm such that the spin configurations with 
potential energies that are closer to the target energy are 
preferentially accepted. In other words, we steer the ini- 
tial state to some target spin configuration such that we 
obtain a state whose potential energy is equal to a specific 
energy. In other respects (i.e., except for the constraint 
of the specified energy), the spin configuration is kept 
random. 

The Hamiltonian equation for Eq. © reads 



1 



cos 



i) 



(2) 



dU 



(3) 
(4) 



where, J is the interaction strength, 9i is the phase angle 
of the spins at site i and the sum is over nearest neigh- 
bor pairs. In the kinetic term, nij = 9i represents the 



where j denotes the nearest neighbors of site i with 
i = 1, ■ • • ,N 2 . For simplicity, we set the rotational in- 
ertia to be equal to unity and put J = 1. Here, we 
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note that, in addition to the conservation of total energy, 
another conserved quantity exists, namely the total ro- 
tational momentum ^ m, = J2i This is due to the 
invariance of the equations under global rotations in the 
variable 9i, that is, under 9i — > 9i + a. 

Equations ([2]) and ([3]) are numerically integrated in 
time using a second order velocity- Verlet algorithm [38| 
with the time integration step of At = 0.01, which con- 
serves the total energy to within a ratio of 10~ 4 up to 
the maximal time integration steps of around 10 6 . Peri- 
odic boundary conditions on both lattice directions are 
employed. The simulations are carried out on square 
lattices of dimensions up to N x N = 2000 x 2000 with 
sample averages of over 40 to 150 different initial configu- 
rations. A parallel computation is employed via a domain 
decomposition method with high parallel efficiency. 
The main quantities of interest are as follows: 
(i) Equal-time spatial correlation function of the order 
parameter, 



■(*) 



(5) 



where < • • • > denotes an average over random initial 
configurations. 

(ii) The total number of vortices and antivortices 
Nv(t) at time t, and 

(iii) Nonequilibrium spin autocorrelation function 



4(^1^0^(0)-^)) 



(6) 



III. SIMULATION RESULTS AND 
DISCUSSIONS 



A. Nonequilibrium relaxation (NER) from fully 
ordered state 

As was mentioned in the introduction, the equal-time 
spatial correlation function in the coarsening process is 
expected to obey a critical dynamic scaling (Eq. (TTJ). 
One of the parameters to be determined is the equilib- 
rium correlation exponent r](E). We may try to fit the 
equal-time spatial correlation function for the coarsening 
dynamics with rj(E) as a free parameter. In this case, 
various values of rj(E) should be chosen on a trial and 
error basis, and the best dynamic scaling collapse must 
be chosen among these many possibilities. This is not 
a very accurate or reliable procedure let alone being te- 
dious. 

Instead, we have a better way to estimate the values of 
rj(E) independently. This is the so-called NER method 



where, in the usual case, the system is instatntaneously 
brought from a perfectly ordered state to a finite tem- 
perature T. This target temperature is usually chosen as 
the critical temperature corresponding to a second order 
phase transition. Here, in our case of Hamiltonian dy- 
namics, the energy per site E is conserved. It is easy to 
see that an analog of NER method here can be achieved 
as follows. Similar to the case of usual lattice spin system 
under thermal noise, here, one starts with a perfectly or- 
dered initial state. But we have to remember that the 
per-site energy of the system has to be conserved (equal 
to E). Since the perfectly ordered initial state has zero 
potential energy (from our definition of the Hamiltonian 
in Eq. ([5])), the rotational velocities 9\ in the initial con- 
figuration should be chosen such that the initial per-site 
kinetic energy be equal to E. Except for this constraint 
of energy conservation, the initial velocities are chosen 
randomly. Beginning with this initial state of perfect or- 
der but with finite kinetic energy (through random rota- 
tional velocities) , when the initial kinetic energy per site 
is low enough, the system will evolve toward a steady 
state corresponding to a critical equilibrium state at a 
finite temperature (correspondng to temperatures below 
the BKT transition). 

By analyzing the approach to the steady state in terms 
of the equal-time spatial correlation function of the order 
parameter, we can compute more accurately the equilib- 
rium correlation exponent r)(E). The NER equal-time 
spatial correlations for the case of E = 0.4 are shown in 
Fig. la. As time elapses, spatial extent of critical power 
law relaxation increases with a clear power exponent. In 
order to extract the r\ exponent more rigorously, we an- 
alyze the data with a critical dynamic scaling 

C nr (r, t) = r-^g nr (r/L nr (t)), (7) 

with the spatial correlation exponent rj{E) and a grow- 
ing correlation length scale L nr (t) with L nr {t) ~ t l l Znr . 
A good collapse of the rescaled correlation functions is 
shown in Fig. lb. Interestingly, we find that the correla- 
tion length scale L nr (t) exhibits a linear growth in time 
as L nr {t) r~j t i.e., l/z nr ~ 1 which is quite different from 
the case of nonconserved dissipative dynamics where the 
diffusive exponent 1/2 is observed [39]. This may be at- 
tributed to the effect of propagating spin wave mode in 
the linearized equation of motion. Indeed, in the spin- 
wave approximation, we find that the equal-time spatial 
correlation C nr (r, t) can be expressed in closed form (see 
the Appendix A) which exhibits critical dynamic scaling 
and linear growth of the length scale L nr (t) with 



a 



l nr:SW (r, t) ~ r-^(E) gnr Sw{r/2t)) (g) 



where f]sw(E) — E/2it and 



gnr-.sw(u) = vP sw 



Vsw(E) 



l + VT 



Vsw(E)/2 



1-vT 



L nsw(E) ( u > ^ 



(u < 1), 



(9) 
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where SW in the subscript refers to "spin-wave approx- 
imation". The behavior of C nr (r,t) vs. r at the latest 
time t — 320 in simulation for several cases of the en- 
ergy E are shown in Fig. lc with the i](E) exponents and 
the spin- wave approximation values rjsw(E) in the inset. 
We can see an approximate agreement between the sim- 
ulation results of r](E) and rjsw(E) especially at low E 
which is consistent with the spin-wave approximation. If 
we consider the case of E = 0.4, the simulation results of 
rj(E = 0.4) = 0.071 can be compared reasonably with the 
approximate analytic value of rjsw(E) = E/2it ~ 0.064 
at spin- wave approximation. Also, with the same approx- 
imation, we could easily derive the formula for the auto- 
correlation A nr . S w{t) as A nr:S w(t) ~ f-»7sw(£0/2 ( see 
APPENDIX A). 



B. Critical coarsening dynamics 

Now we turn to the critical coarsening dynamics evolv- 
ing from random configurations. Figure 2a shows the 
equal-time spatial correlations for E = 0.4 for different 
time instants. We attempt a critical dynamic scaling 



C(r,t)=r-^g(r/L(t)), 



(10) 



where we use the spatial correlation exponent rj(E) ob- 
tained from the NER analysis shown above. In our case, 
for a given time instant t, we determined L{t) in such 
a way that r^C(r, t)\ {r=L{t]) = g(r / L(t))\ [r=L[t)] = 
g(l) = 0.4. Shown in Fig. 2b is such an attempt at crit- 
ical dynamic scaling where we find that the equal-time 
spatial correlation functions exhibit a reasonable critical 
dynamic scaling, at least, in the late-time regime. 

In the late-time regime, we find that the length scale 
L(t) grows with an effective exponent 1/z that is appar- 
ently larger than the diffusive exponent 0.5. For example, 
in the case of E = 0.4, we obtain 1/z ~ 0.553 ± 0.005, 
which is a little larger than 1/2 (Fig. 2c). This is con- 
sistent with the previous work |17| . In contrast, in the 
conventional purely diffusive case, the effective growth 
exponents obtained numerically are invariably smaller 
than 1/2 due to a negative logarithmic correction. Simi- 
lar enhancement of growth was observed recently in the 
coarsening dynamics of the (soft spin) model E as well 
as the MCXY model. There, it was interpreted as the 
result of a logarithmically divergent mobility of vortices 
(and antivortices) in terms of a simplified kinetic model of 
a single pair of vortex and antivortex. Here we find that 
the same model can explain the coarsening process only 
with a modification of the initial condition that accounts 
for the finite length scale Lq at time t = due to the spe- 
cial nature of Hamiltonian dynamics and corresponding 
energy conservation in the present system. Similar phe- 
nomenon of logarithmically divergent mobility was also 
found experimentally in the quasi-two-dimensional diffu- 



sion of colloids 40] or diffusion of protein molecules on 
the membranes 4l| under hydrodynamic effect, and also 
in numerical simulations of the vortex diffusion in the 



anisotropic Heisenberg system in two dimensions with 
spin precession |42h44| ■ 

We reiterate briefly on the single vortex-antivortex an- 
nihilation model. We can see that, in order for the system 
to reach a state with average length scale of R, vortex 
pairs of sizes on the order of R must be already annihi- 
lated. Ignoring the very complicated many-body process 
involving many vortex pairs during this process, we sim- 
plify the whole coarsening process (corresponding to the 
growth of the length scale up to R) by the annihilation 
of single vortex-antivortex pair of size R with a suitably 
defined interaction potential. 

We assume that the vortex acts like a small particle 
with finite 'mass' moving under the influence of an ex- 
ternal force with a mobility that depends logarithmically 
on the length scale of the system. Then the distance R 
between a vortex and an anti-vortex would be described 
by the following equation of motion 



m{R) — 



1 dR 



F(R) = - 5 . 



(11) 



Here, the length-dependent effective mass of a vortex 
is denoted as m(R) and the length-dependent mobil- 
ity of a vortex as /i(i?). We also assume that the 
vortex-antivortex pair is interacting via a Coulombic 
force F(R) = —k/R in two dimensions. We will set the 
proportionality constant k as k = 1. As for the functional 
form of the m(R) and (J-(R), we use the same forms as 
that employed in [5J| where logarithmic dependences of 
the vortex mobility /i(i?) and m(R) are used [42|-[44j as 
follows 

m{R) = to + n%i ln(i?/r ), n(R) = /i + Mi m (-RAo) 

(12) 

where mo, mi, /xo, [i\ are constants and ro denotes a 
shortest cutoff length scale in the system (corresponding 
to the vortex core size). Note that too corresponds to 
the effective mass at the shortest cutoff length scale and 
similarly for /io for the vortex mobility. 

The growth law is obtained by calculating (through in- 
tegration of the above model dynamic equation) the time 
r it takes for a vortex-antivortex pair of size R — R* to 
reach a state with R = Rq where Rq denotes the ini- 
tial (coarsening) length scale of the system (the average 
separation between vortices or antivortices). We set the 
value of dR/dt at R = R* to be zero. We numerically 
solve the model equation of vortex-antivortex pair dy- 
namics, Eqs. (fTTj) and (TT2"j) . to obtain the time r when 
the size of the vortex-antivortex pair R becomes equal to 
Ro- Note that Rq is larger than the lower cutoff length 
scale of ro corresponding to the unit lattice spacing. By 
plotting the resulting relation R* and r we get the growth 
law of the coarsening dynamics. 

If we assume that the inertial effect is negligible, then, 
due to the logarithmic divergence of the vortex mobility, 
we can analytically get the dominant asymptotic growth 
law as 



L(t) ~ (tint) 



1/2 



(13) 
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However, we note that, in the limit of t — >■ 0, the do- 
main length scale L obtained from our simulations of 
our Hamiltonian XY model approaches a finite value 
which is incompatible with the above form (which gives 
Lit = 0) = 0). In order to reconcile with this limiting 
conditions, we should take a more general ansatz equa- 
tion for the time dependence of the length scale that ex- 
hibits a finite initial value for the length scale and still 
satisfies the dynamic equation of motion asymptotically. 
One plausible candidate is the following form which is 
obtained by a simple shift of the time variable t as 

{(t + t )\n(t + t )) 1/2 . (14) 

We might also take another form by adding a constant 
L to Eq. (TT31) as L(t) ~ L + a(Hni) 1 / 2 . But here we 
chose the above method of shifting the time because it 
appears the most natural, being analytic at t = 0. We 
find that L(t) vs. t can be fitted by Eq. (fT4")) in the late 
time regime (Fig. 2c). 

In addition, shown in Fig. 3a is the relaxation of the 
vortex number density p v {t) for E — 0.4. It exhibits 
a power law relaxation with p v (t) ~ t -1,0 ' 6 . If we as- 
sume that the vortices are uniformly (and randomly) dis- 
tributed in two dimensional space, then using a naive 
argument, we would expect that p v (t) ~ L~ 2 {t). Inter- 
estingly, however, we find that the following nontrivial 
scaling relation holds in the late time region (Fig. 3b) 

p v {t) ~ L- X (t), a; =1.93 (15) 

This feature is also very similar to the behavior of MCXY 
model and also that of the model E Q . We note that the 
exponent x showed little systematic variation at different 
values of E except for some statistical numerical fluctu- 
ations (data not shown). 

The nonequilibrium spin autocorrelation function A(t) 
is expected to be related to the growing length scale Lit) 
through a new non-equilibrium exponent A as 

A(t) ~ L-\t). (16) 

We could extract the value of A by plotting A(t) versus 
Lit) as shown in the Fig. 4, where we can see that, in the 
long time limit the value of A approaches A ~ 2.21 for the 
case of E = 0.4. This value is rather close to the value in 
Ref. Q (where A ~ 1.99 was observed) but much larger 
than the value of A ~ 1.17 for the case of non-conserved 
0(2) model with no reversible mode coupling (25l. [32Tj35| . 
It was argued in Ref. Q that the higher mobility of the 
vortices in the present model causes a faster memory loss 
of the initial configuration and hence larger value of the A 
exponent. It would be interesting if one can find some an- 
alytic argument for the value of A = 2.21. It appears that, 
due to the existence of propagating spin wave modes, the 
spin configuration loses its memory more quickly than in 
the case of purely dissipative dynamics. We also found 
that there exists a very weak monotonic increase of the A 
exponent with the increase of E [391 ] . which probably is 



related to the increasing spatial critical exponent r](E) at 
finite temperature. This trend is also consistent with the 
small difference between the value of A exponent in the 
present work and that of MCXY or model E in Ref. Q , 
since, in the latter case, they are dealing with a zero tem- 
perature coarsening. Figure 5 shows a fitting of the full 
simulation results by our vortex-antivortex annihilation 
model via integration of Eq. pT|) with a suitable choice 
of the parameters. 

We note again that, in our Hamiltonian dynamics sim- 
ulations, the initial states are not completely random 
states due to the finite energy constraint, Therefore, the 
initial states already have some finite correlation length 
scale, i.e., separation Rq between vortices. These length 
scales are of the order of 3 to 5 lattice constants. There- 
fore, we expect that a finite time scale To (correspond- 
ing to the time scale for a vortex to travel the distance 
Rq) exists after which a scaling region emerges. It is in- 
teresting to note that the effect of the vortex inertia is 
observed to be not very significant at least for the do- 
main growth law. This might be due to the finite initial 
length scale Lq. This should be contrasted to the case of 
the domain growth in the 0(2) Ginzburg Landau model 
with reversible mode coupling. The typical length scale 
L v (t — 0) between vortices at initial time can be de- 
rived from the number density of vortices, which can be 
compared with the corresponding length scale Lq derived 
from the equal-time spatial correlation functions. We ex- 
pect that these two length scales are approximately pro- 
portional to each other, which is indeed what we found 
from our data [39|. Note that Lq was obtained by the 
arbitrary choice of 0.4 for the value of the cut in the 
dynamic scaling collapse of the equal-time spatial cor- 
relation functions. Since L v (t = 0) represents the ini- 
tial average separation between vortices and antivortices 
at which the value of the rescaled correlation function 
would most probably be smaller than 0.4, we expect that 
L v (t = 0) is larger than Lq. We could also confirm this 
from the data (data not shown here). 



IV. SUMMARY 

We have studied the coarsening dynamics of the Hamil- 
tonian XY model on two-dimensional square lattice. An 
initial state that is specially tuned to have a given po- 
tential energy (otherwise random) but with zero kinetic 
energy develops into a late-time coarsening state, where 
the potential energy slowly decays as a power law with 
a compensating increase in the kinetic energy. In order 
to compute accurately the equilibrium correlation expo- 
nents r)(E) for the order parameter, we employed an ana- 
log of the so-called NER method where the system be- 
gins with a perfect order but with finite kinetic energy. 
An interesting scaling behavior emerged here with lin- 
ear growth of the nonequilibrium correlation length. In 
the coarsening dynamics, the growth of the characteristic 
length scale can be fitted as Lit) ~ ((t + to) ln(t + to)) 1 ^ 2 
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with a finite offset time to which is due to the existence 
of a finite separation between vortices and antivortices at 
the initial time. The growth of the length scale can be in- 
terpreted as the vortex mobility diverging logarithmically 
in the size of the vortex-antivortex pair. A simple phe- 
nomcnological dynamic model of vortex-antivortex anni- 
hilation can fit the growth of the length scale, with a suit- 
able set of parameters for the effective inertia and mobil- 
ity of vortices that are dependent on the growing length 
scale. The nonequilibrium spin autocorrelation function 
A(t) and the L(t) are related through the nonequilib- 
rium exponent A by A{t) ~ L~ x {t) with an exponent of 
A ~ 2.21 (at E = 0.4) which is distinctly different from 
the case of ordinary 0(2) or XY model with dissipative 
dynamics. 
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Appendix A: Critical dynamic scaling of the 
correlation functions in spin-wave approximation 

Here we analytically show that the critical dynamic 
scaling Eq. (7) holds in the spin-wave approximation. 
We consider the equation of motion for the spin-angle 
variables 8i . For convenience of notation, we denote here 
9i as 9(r, t) at position r (instead of i) and time t. In the 
spin-wave approximation with continuum notation, the 
Hamiltonian equation of motion reads as follows 



d 2 9(r,t) 
dt 2 



V 2 d(r,t)=0 



(Al) 



which is the same as a wave equation. Here we simply 
put the wave velocity c = 1. We want to calculate the 
equal-time spatial correlation of the spins for the case 
of perfectly ordered initial spin configuration (but with 
finite angular velocities) which is defined as 



C nr {r, t) = (exp [i0(r, t) - i0(0, t)}) 



(A2) 



and the nonequilibrium spin autocorrelation function 

A nr (t) = (exp [i6{r, t) - i0(r, 0)]) . (A3) 

General solutions for the wave equation can be written 
as 

(A4) 



where J k = fd 2 k/(2ir) 2 and u> = |k| 
initial condition of 9(r,t = 0) = 0, we get f£ = — = 
/k/2i with 



k. With the 



6(r,t) = / /kcxp(ik • r)sin(wi). (A5) 
Jk 

We assume that the initial angular velocity at each site 
has a Gaussian distribution with zero mean and stan- 
dard deviation proportional to the square root of the av- 
erage energy per site. In addition, it is also assumed that 
the angular velocities exhibit delta-function correlations 
between different spatial points. If the average energy 

per site is given by E, then we have (o{r, O)0(r', 0)^ = 

2E5(r — r'), which, in the wave- vector space, can be writ- 
ten as 



(AW = ^(k + k'). 



(A6) 



Combining above relations and also using the gaussian 
nature, we can obtain the following 



a 



nr:SW 



(r,t) 



exp 



= exp 



-±((0(r,t)-0(O,t)) 2 



E 

^2 



d 2 k -^sin 2 (ik-r) sin 2 (fct) 
(A7) ' 



In order to compute the integral inside the exponential, 
we put 

J(f, t) = J d 2 k i sin 2 (Ik • r) sin 2 (/rf) (A8) 

Taking a derivative of J with respect to t and also adding 
prescription for the proper convergence in the (positive) 
large k limit, we arrive at 



dJ(r,t) 



= lim e ^ +0 \ Jq 71 def) / °° dk sin(2/rf) 



dt ' ' (A9) 

x ( 1 — cos (fc r cos 4>) J exp(— ek). 

This integral can now be easily done with the following 
result 



dJ(r,t) 
dt 



Mi-p-W^-v) { L t<1) , 



It 
7T 

2t 



(A10) 

The above equation can now be integrated with respect to 
t with proper boundary conditions, leading to the equal- 
time spatial correlation function as 



C nr:S w(r,t) ~ r-^( E ) 9nr:SW (r/2t) 



(AH) 
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where T]sw(E) = E/2ir and 



9nr:Sw(u) = U^( E ) ^ ^ _ ^ j (« < 1), 

= u nsw(E) ( u > X y 

(A12) 

In the limit of u = r/2t << 1, we get C nr: sw(r, t) <~ 
f-isirfE)^ _ ^ SH/ (i^) M 2 /8) and also in the other regime 
of u ee r/2t > 1, we get C nr:SW (r,t) ~ r''™^ in- 
dependent of the value of r. This agrees qualitatively 



with the simulation results including the plateau behav- 
ior in the region of u = r/2t > 1. The value of rjsw{E) 
at spin- wave approximation for E = 0.4 is given by 
E/2ir ~ 0.064 which is in approximate agreement with 
the simulation result of rj(E = 0.4) = 0.071 (see Fig. lc). 
Also we find for the autocorrelation A nr: sw{t) 

A nr:SW (t) ~ t-"^^)/ 2 (A13) 

which again shows approximate agreement with simula- 
tions results (data not shown). 



[1] J. D. Gunton, M. San Miguel, and P. S. Sahni, in Phase 
Transitions and Critical Phenomena, edited by C. Domb 
and J. L. Lebowitz (Academic, New York, 1983), Vol 8; 
H. Furukawa, Adv. Phys. 34, 703 (1985); K. Binder, Rep. 
Prog. Theor. Phys. 50, 783 (1987). 

A. J. Bray, Adv. Phys. 43, 357 (1994). 
S. Ma and G. F. Mazenko, Phys. Rev. B 11, 4077 (1975). 
D. R. Nelson and D. S. Fisher, Phys. Rev. B 16, 4945 
(1977). 

K. Nam, B. Kim, and S. J. Lee, J. Stat. Mech. Theo. 
Exp., 2011, P03013 (2011). 

P. C. Hohenberg and B. I. Halperin, Rev. Mod. Phys. 49, 
435 (1977). 

R. Folk and G. Moser, J. Phys. A: Math. Gen. 39, R207 
(2006). 

L. Casetti, M. Pettini, and E. G. D. Cohen, Phys. Rep. 
337 237 (2000). 

L. Caiani, L. Casetti, C. Clementi, G. Pettini, M. Pettini, 
and R. Gatto, Phys. Rev. E 57, 3886 (1998). 
L. Caiani, L. Casetti, and M. Pettini, J. Phys. A 31, 
3357 (1998). 

X. Lconcini, A. D. Verga, and S. Ruffo, Phys. Rev. E 
57, 6377 (1998). 

S. Lepri, and S. Ruffo, Europhys. Lett. 55, 512 (2001). 
M. Cerruti-Sola, C. Clementi, and M. Pettini, Phys. 
Rev. E 61, 5171 (2000). 

V. Latora, A. Rapisarda, and S. Ruffo, Physica D 131, 
38 (1999). 

B. Zheng, Phys. Rev. E 61, 153 (2000). 
J. Kockelkoren and H. Chate, Phys. Rev. E 65, 058101 
(2002). 

K. Koo, W. Baek, B. Kim, and S. J. Lee, J. Korean Phys. 
Soc. 49, 1977 (2006). 

V. L. Berezinskii, Zh. Eksp. Teor. Fiz. 59, 907 (1970) 
[Sov. Phys. JETP 32, 493 (1971)]; J. M. Kosterlitz, and 
D. J. Thouless, J. Phys. C 6, 1181 (1973); J. M. Koster- 
litz, ibid. 7, 1046 (1974). 

R. Loft and T. A. DeGrand, Phys. Rev. B 35, 8528 
(1987); H. Toyoki and K. Honda, Prog. Theor. Phys. 
78, 237 (1987). 

H. Toyoki, Phys. Rev. A 42, 911 (1990). 
M. Mondello and N. Goldenfeld, Phys. Rev. A 42, 5865 
(1990); Phys. Rev. E 47, 2384 (1993). 

A. J. Bray and K. Humayun, J. Phys. A 23, 5897 (1990); 
S. Puri and C. Roland, Phys. Lett. A 151, 500 (1990); 
H. Toyoki, Phys. Rev. A 42, 911 (1990). 

B. Yurke, A. N. Pargellis, T. Kovacs, and D. A. Huse, 



Phys. Rev. E 47, 1525 (1993). 
[24] R. E. Blundell and A. J. Bray, Phys. Rev. E 49, 4925 
(1994). 

[25] J-R. Lee, S. J. Lee, and B. Kim, ibid. 52, 1550 (1995). 
[26] F. Rojas and A. D. Rutenberg, Phys. Rev. E 60, 212 
(1999). 

[27] A. J. Bray, A. J. Briant, and D. K. Jervis, Phys. Rev. 

Lett. 84, 1503 (2000); 
[28] A. J. Bray, Phys. Rev. E 62, 103 (2000). 
[29] H. P. Ying, B. Zheng, Y. Yu, and S. Trimper, Phys. 

Rev. E 63, 035101 (2001). 
[30] A. A. Ovchinnikov and Ya. B. Zeldovich, Chem. Phys. 

28, 214 (1978); D. Toussaint and F. Wilczek, J. Chem. 

Phys. 78, 2642 (1983); K. Kang and S. Redner, Phys. 

Rev. A 30, 2833 (1984). 
[31] G. S. Jeon, S. J. Lee, and M. Y. Choi, Phys. Rev. B 67, 

014501 (2003). 

[32] T. J. Newman and A. J. Bray, J. Phys. A 23, L279 (1990); 

23, 4491 (1990). 
[33] A. J. Bray and K. Humayun, J. Phys. A 23, 5897 (1990); 
[34] T. J. Newman, A. J. Bray, M. A. Moore, Phys. Rev. B 

42, 4514 (1990). 
[35] F. Liu and G. F. Mazenko, Phys. Rev. B 45, 6989 

(1992); 46, 5963 (1992). 
[36] J. Das and M. Rao, Phys. Rev. E 57, 5069 (1998). 
[37] J. Das and M. Rao, Phys. Rev. E 62, 1601 (2000). 
[38] H. Yoshida, Phys. Lett. A 150, 262 (1990). 
[39] K. Nam, S. J. Lee, and B. Kim, unpublished. 
[40] J. Sane, J. T. Padding, A. A. Louis, Phys. Rev. E 79, 

051402 (2009). 

[41] P. G. Saffman and M. Delbriick, Proc. Nat. Acad. Sci. 

72, 3111 (1975). 
[42] F. G. Mertens and A. R. Bishop, Dynamics of Vortices in 

Two-Dimensional Magnets, editted by P. L. Christiansen 

and M. P. Sorensen, Nonlinear Science at the Dawn of 

the 21st Century, (Springer, Berlin, 1999). 
[43] T. Kamppeter, F. G. Mertens, E. Moro, A. Sanchez, and 

A. R. Bishop, Phys. Rev. B 59, 11349 (1999). 
[44] T. Kamppeter, F. G. Mertens, A. Sanchez, A. R. Bishop, 

F. Dommguez-Adame, and N. Gr^nbech- Jensen, Eur. 

Phys. J. B 7, 607 (1999). 



9 



1 






FIG. 1: (a) The equal-time spatial correlation function of the 
XY spin order parameter with NER dynamics for E — 0.4 
at times t — 10, 20, 40, 80, 160, 320, where the system size 
is 1800 x 1800. (b) The scaling collapse of the data in (a) 
with the appropriate scaling length L„ r (t). The inset shows 
the length scale L nr (t) vs. t which exhibits a linear growth 
L nr (t) ~ t 1/z with 1/z ~ 1 and rj ~ 0.071. (c) The equal-time 
spatial correlations C nT (r,i) at late time t = 320 for different 
energies. The inset shows the corresponding r/(E) (circles) 
and the spin- wave approximation values r)sw(E) (solid line) 
versus E. 





FIG. 2: (a) The equal-time spatial correlation function of 
the XY spin order parameter with coarsening dynamics for 
E — 0.4 at various times from t = 10 up to 5120 with the 
system size is 2000 x 2000. (b) The scaling collapse of the data 
in (a) with the appropriate scaling length L(t) for t — 160, 
320, 640, 1280 and 2560. The inset shows the length scale 
L(t) vs. t which exhibits an effective superdiffusive power law 
growth of L(t) ~ t 1/z with 1/z ~ 0.553. (c) L 2 (t)/(t + t ) vs. 
ln(i + fo) which indicates a positive logarithmic correction for 
the growing length scale. 




FIG. 3: Relaxation of (a) the vortex number density p v (t) 
at energy E = 0.4 and system size 2000 x 2000 which shows 
a power law decay with p v (t) ~ £ _1 076 ( m the time interval 
between t = 10 and 1000). (b) The vortex number density 
pv{t) vs. L(t) exhibits a power law decay p v (t) ~ L~ x (t) with 
x ~ 1.93. 




FIG. 4: Spin autocorrelation function vs. L(t) at energy 
E — 0.4 indicating a simple power law relationship between 
the two with A(t) ~ L(t)~ 2 21 . 
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FIG. 5: The growth of L(t) fitted by a simple vortex- 
antivortex annihilation model at energy E = 0.4 which shows 
a good agreement. 



